## This program generates Figure 10 from "Financial Constraints, Sectoral Heterogeneity,
##  and the Cyclicality of Investment" by Cooper Howes

##### Step 0: Preliminaries #####

## Clear variables and data
rm(list=ls())

## Load packages
library(tidyverse)

##### Step 1: Load data #####

## Baseline IRFs 
## Created by dynare file "baseline_model.mod"
irfdata <- t(read.csv("irfs/shock_irfs_baseline.csv",header=F)[1:9,])
colnames(irfdata) <- irfdata[1,]
rownames(irfdata) <- NULL
irfdata <- irfdata[-1,]
irfdata <- apply(irfdata, 2, as.numeric)
irfdata <- irfdata[,c("k_d_e_m","k_n_e_m","k_e_m")]

## IRFs for model with fully flexible prices in both sectors 
## Created by Dynare file "flexprice_model.mod"
irfdata_flexprice <- t(read.csv("irfs/shock_irfs_flexprice.csv",header=F))
colnames(irfdata_flexprice) <- irfdata_flexprice[1,]
rownames(irfdata_flexprice) <- NULL
irfdata_flexprice <- irfdata_flexprice[-1,]
irfdata_flexprice <- apply(irfdata_flexprice, 2, as.numeric)

## IRFs for model without adjustment costs
## Created by Dynare file "nofrictions_model.mod"
irfdata_frictionless <- t(read.csv("irfs/shock_irfs_frictionless.csv",header=F)[1:9,])
colnames(irfdata_frictionless) <- irfdata_frictionless[1,]
rownames(irfdata_frictionless) <- NULL
irfdata_frictionless <- irfdata_frictionless[-1,]
irfdata_frictionless <- apply(irfdata_frictionless, 2, as.numeric)

## IRFs for model without firm financial constraints
## Created by Dynare file "nofirmconstraints_model.mod"
irfdata_nofirmconstraints <- t(read.csv("irfs/shock_irfs_nofirmconstraints.csv",header=F)[1:9,])
colnames(irfdata_nofirmconstraints) <- irfdata_nofirmconstraints[1,]
rownames(irfdata_nofirmconstraints) <- NULL
irfdata_nofirmconstraints <- irfdata_nofirmconstraints[-1,]
irfdata_nofirmconstraints <- apply(irfdata_nofirmconstraints, 2, as.numeric)




##### Step 2: Plot alternative model specifications #####

## Create x-axis labels
qtrs <- 0:16

## Create new window
windows(width=40,height=25)
par(oma=c(5.5,2,1,1))
m <- matrix(1:4,nrow=2,byrow=T)
layout(mat = m)

marpar <- c(2,3,3,1)

## Specify y-axis range
range_cap2 <- c(-0.02,0.06)

## Capital plots
# Actual (Note: uses same range as data for effect)

# My model
par(mar = marpar)
plot(x=qtrs,y=irfdata[,1],col="blue",type="b",pch=1,lwd=4,ylim=range_cap2,cex.axis=1,
     ylab="",cex.lab=1,xlab="",axes=F,frame=T)
axis(1,cex.axis=1.5)
axis(2,cex.axis=1.5,at=pretty(range_cap2),lab=paste0(pretty(range_cap2)*100, "%"),las=TRUE)
lines(x=qtrs,y=irfdata[,2],type="b",pch=2,col="red",lwd=4)
lines(x=qtrs,y=irfdata[,3],type="b",pch=3,col="orange",lwd=4)
abline(a=0,b=0,col="black",lwd=1)
title(main="Baseline model", col.main="blue", font.main=2,cex.main=1.5,line=1)

## Flexible prices
par(mar = marpar)
plot(x=qtrs,y=irfdata_flexprice[,1],col="blue",type="b",pch=1,lwd=4,ylim=range_cap2,cex.axis=1,
     ylab="",cex.lab=1,xlab="",axes=F,frame=T)
axis(1,cex.axis=1.5)
axis(2,cex.axis=1.5,at=pretty(range_cap2),lab=paste0(pretty(range_cap2)*100, "%"),las=TRUE)
lines(x=qtrs,y=irfdata_flexprice[,2],type="b",pch=2,col="red",lwd=4)
lines(x=qtrs,y=irfdata_flexprice[,3],type="b",pch=3,col="orange",lwd=4)
abline(a=0,b=0,col="black",lwd=1)
title(main="Flexible prices", col.main="blue", font.main=2,cex.main=1.5,line=1)

## No bells/whistles
par(mar = marpar)
plot(x=qtrs,y=irfdata_frictionless[,1],col="blue",type="b",pch=1,lwd=4,ylim=range_cap2,cex.axis=1,
     ylab="",cex.lab=1,xlab="",axes=F,frame=T)
axis(1,cex.axis=1.5)
axis(2,cex.axis=1.5,at=pretty(range_cap2),lab=paste0(pretty(range_cap2)*100, "%"),las=TRUE)
lines(x=qtrs,y=irfdata_frictionless[,2],type="b",pch=2,col="red",lwd=4)
lines(x=qtrs,y=irfdata_frictionless[,3],type="b",pch=3,col="orange",lwd=4)
abline(a=0,b=0,col="black",lwd=1)
title(main="No adjustment frictions", col.main="blue", font.main=2,cex.main=1.5,line=1)


# No firm financial constraints
par(mar = marpar)
plot(x=qtrs,y=irfdata_nofirmconstraints[,1],col="blue",type="b",pch=1,lwd=4,ylim=range_cap2,cex.axis=1,
     ylab="",cex.lab=1,xlab="",axes=F,frame=T)
axis(1,cex.axis=1.5)
axis(2,cex.axis=1.5,at=pretty(range_cap2),lab=paste0(pretty(range_cap2)*100, "%"),las=TRUE)
lines(x=qtrs,y=irfdata_nofirmconstraints[,2],type="b",pch=2,col="red",lwd=4)
lines(x=qtrs,y=irfdata_nofirmconstraints[,3],type="b",pch=3,col="orange",lwd=4)
abline(a=0,b=0,col="black",lwd=1)
title(main="No firm financial constraints", col.main="blue", font.main=2,cex.main=1.5,line=1)


par(fig = c(0, 1, 0, 1), oma = c(0, 0, 0, 0), mar = c(0, 0, 0, 0), new = TRUE)
plot(0, 0, type = "n", bty = "n", xaxt = "n", yaxt = "n")
legend(x="bottom",inset=0,legend=c("Durable","Nondurable","Total"),
       lwd=2,lty=1,pch=c(1,2,3),col=c("blue","red","orange"),cex=1.5,horiz=TRUE)





